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Fast  and  accurate  techniques  for  determining  electron  fluxes  and  energy  deposition 
profile  in  complex  multidimensional  microcircuits  have  been  developed.  The  multi¬ 
dimensional  Spencer-Levis  Transport  equation  is  solved  deterministically  using 
several  special  numerical  techniques. 

We  begin  with  the  derivation  of  the  Spencer-Lewis  equation.  Then  we  discuss  SMART 
scattering  theory  that  enables  us  to  replace  the  highly  anisotropic  electron  scatter¬ 
ing  kemal  by  one  that  is  more  amenable  to  numerical  treatment.  Next,  we  describe 
several  transport  and  diffusion  theory  solution  algorithms.  Finally,  we  present 
new  analytical  benchmarking  methods  that  will  prove  useful  in  generating  more  com¬ 
prehensive  benchmarks. 
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MULTIDIMENSIONAL  DETERMINISnC  ELECTRON  TRANSPORT  CALCULATIONS 


W.  L.  Fnippoae,  B.  D.  Ganapol,  S.  P.  Monahan 


Introduction 

This  report  summarizes  research  carried  out  for  the  US  Air  Force  Rome  Air  Development 
Center  under  contracts  F  30602-81-C-0185  (P.  R.  No.  S-6-7548)  and  E-21-669-S6  (P.  R.  No.  S- 
7-7322).  The  object  of  this  work  was  to  develop  fast  and  accurate  techniques  for  determining 
electron  fluxes  and  energy  deposition  profiles  in  complex  multidimensional  microcircuits,  for  use  in 
radiation  hardness  studies.  Our  approach  was  to  use  deterministic  solutions  of  the  multidimensional 
Spencer-Lewis  electron  transport  equation.  To  do  this  it  was  necessary  to  develop  several  special 
numerical  techniques  which  we  describe  below. 

We  begin  with  the  deriviation  of  the  Spencer-Lewis  equation.  Then  we  discuss  SMART 
scattering  theory  that  enables  us  to  replace  the  highly  anisotropic  electron  scattering  kernel  by  one 
that  is  more  amenable  to  numerical  treatment  Next,  we  describe  several  transport  and  diffusion 
theory  solution  algorithms.  Finally  in  Section  IV  we  present  several  new  analytical  benchmarking 
methods  that  will  prove  useful  in  generating  more  comprehensive  benchmarks. 

I.  Deriviation  of  the  Spencer-Lewis  Equation 

When  continuous  slowing  down  theory  is  valid,  electron  transport  in  homogeneous  media  is 
governed  by  the  Spencer-Lewis  equation, 

[^  +  +  o(s)j  *(r,s,n)  -  Idir  o(s^-*n)  +  QC?4,n)  (D 


where 

4(f  ,s,n)  «  electron  density  (electrons  per  unit  volume  per  unit  solid  angle)  at  position  r,  direction 
n  and  path  length  s. 
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o(s)  -  total  (scattering)  cross  section,  and 

>  electron  source  density  (electrons  per  unit  volume  per  unit  path  length  per  unit  solid 
angte). 

According  to  the  continuous  slowing  down  approximation  (CSDA)  there  is  a  one-to-one 
cmrrespondence  between  distance  traveled  and  energy  lost,  and  the  path  length  s  plays  the  role  of  an 
energy  variable.  We  prefer  the  Spencer-Lewis  equation  to  the  multigroup  Boltzman  equation, 
because  the  latter  is  not  particulariy  efficient  at  modeling  electron  energy  loss.  To  obtain  good 
results  normally  requires  an  exceedingly  large  number  of  groups^**  (or  other  high  resolution  energy 
discretization  technique^),  while  coarse  path  length  mesh  solutions  of  the  Spencer-Lewis  equation  are 
normally  sufficient,^’^  probably  because  the  CSDA  is  built  into  the  equation.  Equation  (1)  was 
presented  by  Lewisf  in  1950  without  proof,  and  its  first  practical  solutions  for  the  infinite  medium 
case  were  obtained  by  Spencer*  in  1955  using  the  method  of  moments.  We  have  derived^  £q.  (1) 
and  a  time-dependent  generalized  Spencer-Lewis  equation  valid  for  inhomogeneous  media  from  the 
Boltzmann  equation. 


+V*V  +  o(r77)v 


where 

N  «  electron  density  (electrons  per  unit  volume  per  unit  velocity)  at  position T,  velocity'^ 

and  time  t, 

o(f  ,v)  a  total  (scattering)  cross  section 

a  differenthd  scattering  cross  section 

qCry,t)  a  electron  source  density  (electrons  per  unit  volume  per  unit  velocity  per  unit  time) 


-  3  - 


A.  Homogeneous  Media 

To  derive  Eq.  (1)  from  Eq.  (2)  we  assume  that 

1.  Scattering  changes  direction  but  not  energy,  that  is, 

-  «(v-v')  o(v.tr-*0)/v*.  (3) 

2.  Energy  loss  is  completely  determined  by  the  stopping  power  | 

3.  The  medium  is  homogeneous. 

4.  The  fixed  source  is  of  the  form 

q(?.‘^.t)  -  Q(?.V)  S{t  -  T(v)l  ,  (4) 

where  T(v)  is  the  time  it  takes  for  an  electron  to  slow  down  to  speed  v  from  some  reference  speed 
Vg,  such  as  the  largest  speed  in  the  system.  (This  source  form  is  needed  so  that  at  time  t  all 
electrons  will  have  the  same  speed  V(t).) 

It  can  then  be  shown  that 

N(r:^,t)  -  N(FAt)  5(v-V(t)l/v*  ,  (5) 

Id  B 

and  Eq.  (1)  restilts  from  the  substitution  -  ^  ^  in  Eq.  (2)  where 

V  eft  OS 


and  Eq  >  E(Vq).  The  details  of  this  deriviation  are  given  in  Ref.  7. 

If  the  source  restriction  [Eq.  (4)]  is  relaxed,  then  the  one-to-one  correspondence  between  t 
and  V  and  therefore  between  t  and  s  is  lost;  however,  Eq.  (1)  remains  valid  as  is  evident  by  the 
alternative  deriviation  discussed  below. 
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B.  Inhomogeneous  Media 

To  obtain  a  Spencer-Lewis  type  equation  for  inhomogeneous  media  we  assume  that 
1.  Energy  loss  and  angular  deflection  are  unconnected,  that  is  that  is  of  the 

form 


o(r^'-»V)  -  «(E'-E)  On(?,E.nr-4i)  +  5(ir-n)  (7fe(?.F-E)  (7) 

2.  An  electron  loses  exactly  the  quantity  of  energy  AE  with  each  collision  while 
preserving  the  correct  stopping  power  so  that 

-E)  -  ^  I  §  I  «(E'.E-AE)  .  (8) 

Inserting  Eqs.  (7)  and  (8)  into  Eq.  (1),  taking  the  limit  AE  -»  0,  and  making  the 
substitution 


JL 

dE 


1 


dE 

ds 


(?.E) 


d_ 

dS 


(9) 


we  obtain  after  considerable  manipulation  [See  Ref.  7] 


— i—  +  c(?3,n)  +  T(f ,s)  i- 

v(?3)  ^ 


^3l,at) 


k 

1 


dir  <7(?3,nr-.n)  ^3.at)  +  q(?3,n,t) 


(10) 


where 


-  5  - 


3(f  3.0)  -  (7(?3)  - 


If 


(11) 


T(?3)  -  1  + 


fi-V  Eff3) 


(12) 


Equation  (10)  is  a  generalization  of  the  Spencer-Lewis  equation  valid  for  spatially 
varying  stopping  powers.  If  the  medium  is  homogeneous  then 


1  -*  s,  (13) 

3(?3.0)  oii)  (14) 


T(?3)  -  1 


(15) 


and  Eq.  (10)  becomes 


1  i. 

V  at 


» 

m  I 


+  t>v  +  o(s)  +  ^  «r,s,0,t)  -  dtr  (Ks.rr-^n)  ^^.ir.t)  q(f,s,at) 


(16) 


Equation  (16)  is  not  simply  the  time-dependent  form  of  Eq.  (1).  The  variable  $  in  the  later 
equation  is  a  number  density  (electrons/cm^)  while  ^  in  the  former  equation  is  a  flux  (cm 
traveled/cm Vsecond).  Integrating  Eq.  (16)  over  all  time  and  putting 

-oo 

«(r,s,n)«  dt^^At)  (17) 

Jo 


and 
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OtfAft)*  dtq(?4At)  (18) 

Jo 

we  obtain  Eq.  (1).  This  deriviatioa  of  £q.  (1)  does  not  require  the  source  restriction  specified  by 
Eq.  (4). 

An  independent  verification^  of  Eq.  (10)  was  obtained  for  the  one-dimensional  case  using 
particle  balance  considerations. 


n.  Numerical  Transport  Methods 

A.  SMART  Scattering  Theory 

Electrons  interact  via  long-range  electric  forces,  and  therefore  their  cross  sections  tend  to  be 
enormous  and  extremely  forward  peaked  as  compared  to  those  of  neutral  particles.  As  a 
consequence,  electron  transport  is  characterized  by  a  large  number  of  minute  deflections  in  the 
velocities  of  the  electrons.  This  virtually  rules  out  the  possibility  of  direct  modeling  either  through 
Monte  Carlo  or  discrete  ordinates  techniques. 

Most  Monte  Carlo  codes  use  condensed  histories***  because  the  large  value  of  the  scattering 
cross  section  makes  it  impractical  to  follow  each  individual  scattering  event  for  a  sufficient  number 
of  histories. 

The  large  scattering  cross  sections  also  cause  trouble  for  discrete  ordinates  codes  by  making 
the  spectral  radius  (without  acceleration)  for  the  source  iteration  nearly  equal  to  \mity.  Furthermore, 
the  anisotropy  of  the  scattering  kernel  can  require  a  large  number  of  discrete  ordinates  for  an 
adequate  numerical  treatment.  Because  of  these  problems,  discrete  ordinates  solvers  normally  replace 
the  true  transport  model  by  one  that  is  more  easily  simulated.  This  is  done  through  the  use  of 
effective  cross  sections.  These  sections  have  three  important  properties  that  make  them  suitable  for 
numerical  calculations; 

-They  are  much  smaller  than  the  true  cross  sections  they  replace. 
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-They  are  much  les  anisotropic. 

-They  yield  good  results  in  discrete  ordinates  codes. 

With  effective  cross  sections  a  few  large  angle  deflections  are  used  to  model  the  combined  effect  of 
many  small  angle  deflections. 

Effective  scattering  matrices  are  most  often  generated  using  Fokker-Planck^o*!^  methods  or 
the  extended  transport  (delta  function)  correction.^^*^^  The  major  drawback  with  these  techniques  is 
that  they  lead  to  nonpositive  and  therefore  nonphysical  scattering  matrices. 

We  describe  a  special  type  of  effective  scattering  matrix  that  we  refer  to  as  SMART 
i[simolation  of  n^y  jccumulative  Rutherford_mjectories).  The  acronym  is  somewhat  appropriate 
(although  not  in  the  precise  artificial  intelligence  sense)  because  the  matrix  elements  are  defined 
such  that  they  cancel  errors  due  to  angular  discretization.  It  is  this  property  that  enables  electron 
discrete  ordinates  calculations  to  be  performed  with  relatively  few  (typically  12)  discrete  directions. 
It  is  also  possible  to  refine  the  definition  of  these  scattering  matrices  such  that  they  cancel  path 
length  in  addition  to  angular  discretization  errors.  We  refer  to  such  scattering  matrices  as  very 
SMART. 

The  theory  of  SMART  scattering  matrices  is  based  on  the  conjecture  that  a  scattering  kernel 
should  be  independent  of  the  problem  geometry.  In  particular,  a  scattering  matrix  that  performs 
well  in  an  infinite  medium  should  also  perform  well  in  a  finite  medium.  Using  the  Goudsmit- 
Saunderson^*  theory  of  multiple  scattering,  exact  infinite  medium  solutions  are  easily  obtained. 
Working  backward  from  such  results,  we  are  able  to  deduce  a  suitable  scattering  matrix  for  the 
discrete  ordinates  equation.  In  contrast  to  the  Fokker-Planck  kernel  and  the  extended  transport 
corrected  kernel  (which  is  SMART  with  Gauss  quadrature  sets  but  not  positive),  positive  SMART 
scattering  matrices  are  easily  generated. 

Positive  scattering  matrices  have  three  important  advantages; 

-They  guarantee  positive  solutions  when  used  in  conjunction  with  positive  spatial 


differencing  schemes. 
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-They  enable  the  use  of  negative  flux  fixups.  (Such  fixups  are  justified  when  negative 
ai^ular  fluxes  result  from  spatial  differencing  errors,  but  not  when  they  resxilt  from  negative 
scattering  matrix  elements.) 

-They  can  be  used  in  single  collision  Monte  Carlo  codes  that  require  positive  cross  sections. 


There  are  infinitely  many  ways  to  discretize  the  electron  transport  equation,  an  infinite 
subset  of  which  leads  to  the  discrete  ordinates  equations.  Corresponding  to  each  discretization 
scheme  is  a  different  SMART  scattering  matrix.  In  fact,  if  Gauss  quadrature  points  are  used,  then 
one  of  these  schemes  leads  to  the  extended  transport  corrected  scattering  matrix.  Fortunately,  it  is 
easy  to  find  discretization  methods  that  lead  to  discrete  ordinates  equations  and  to  positive  SMART 
scattering  matrices  as  weU. 

1.  The  SMART  Scattering  Matrix 

To  illustrate  the  discretization  possibilities,  we  begin  with  a  weighted  residual  deriviation  of 


the  Sff  equations.  (This  section  is  condensed  from  Ref.  IS.) 

The  Weighted  Residual  Deriviation  of  the  Sm  Equations 


Let  the  flux  in  £q.  (1)  be  approximated  by 


-  Q(?4,n) 


(20) 


-  9  - 


Forcing  the  residual  to  be  orthogonal  to  a  set  of  test  functions,  ot  •  1,2,...>1},  that  satisfy 

the  conditions 


A  ^  ^  .A 


T“{n)*B«  (0)  dn  - 


(n)  dn  -  iT»«n 


we  obtain 


+  Q“(r,s)  ,  m  »  1,2,...,M  , 


where 


^  mm*  ”  ^^vaxa  » 


wjj 


A  ^  A  A  ^  #  ^A 


T“(n)*o(s,n'-n)  b«  (ir)  dn*  dn  , 


r,s)  »  J" 


A  ^  .  A, 


crcr.s)  -  T»(n)*Q(r,s,n)  dn  , 


and  from  Eq.  (22), 


A  ^A  A 


nF»  -  T»(n)*nB“‘ (n)  dn 


« n™  , 


(27) 
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wlMre  is  the  m'th  discrete  direction.  The  are  the  elements  of  the  conventional  scattering 
matrix  S'.  The  matrix  S  defined  by  Eq.  (24)  has  removal  as  well  as  scatter  included  in  its  definition. 
It  is  referred  to  as  the  net  conventional  scattering  matrix,  and  is  normally  approximated  as 

-  o(s>*'-»l^)  ,  (28) 

where 

■  |B®(n)  dn  .  (29) 

Because  of  the  extreme  anisotropy  of  o(sjy-*n),  the  major  contribution  to  the  integral  of  Eq. 
(1)  comes  from  values  o(  VC  -il  Unless  a  very  large  value  of  M  is  used,  this  contribution  is 
normally  missed  with  conventional  scattering  matrices.  For  this  reason  the  error,  defined  as 

<®(r,s)  ■  ^(?4)  -  ^(?^)  .  (30) 

can  be  quite  large.  To  reduce  this  error,  we  replace  S  with  a  SMART  net  scattering  matrix  S, 
which  is  designed  to  give  good  results  when  a  coarse  angular  mesh  is  used.  The  matrix  S  is 
obtained  by  forcing  <®  «  0  for  a  set  of  problems  for  which  analytic  solutions  are  available. 

The  Goudsmit-Saunderson  Matrix 

To  obtain  benchmark  problems  that  are  amenable  to  exact  analytical  solution,  we  consider  an 
infinite  source-free  medium  with  an  initial  (s  «  0)  distribution  ^0,0)  that  is  independent  of  position 
T.  For  this  case,  Eq.  (1)  reduces  to 


(31) 
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We  are  interested  in  equating  the  discrete  ordinates  and  exact  solutions  for  identical  initial 
conditions.  As  seen  from  £q.  (19),  the  initial  condition  for  the  discrete  ordinates  equatkms  most  be 
of  the  form 


^(0.n)  -  ^  ^(0)  B"»(n)  .  (32) 

m«l 


With  an  initial  distribution  of  this  form,  the  exact  0™(s)  are  given,^*  in  matrix  form  by 


7(s)  -  G(0-s)7t(0)  . 


(33) 


where  the  components  of  ^  and  are  ^  and  respectively.  The  elements  of  the  Goudsmit- 
Saunderson  matrix  G  are  given  by 


G(So— ^)i 


'0  ^'xmo 


V  a 


T®(n)*g(s--*s,  (in')  B“'(ir)  dn  dtr  , 


where*^ 


g(So-*s,nfr) 


OO 

y  exp.  -1  l(7o(s')  -  0^(3*)]  ds' 


Pi(nir)  , 


(34) 


(35) 


and  P/  is  the  ^th  Legendre  polynomial.  We  emphasize  here  that  Eq.  (33)  represents  the  exact 
(continuous  angle)  solution  for  any  source  that  can  be  expressed  by  Eq.  (32).  Although  spatial 
effects  have  been  eliminated,  the  effects  of  multiple  scattering  are  present  and  therefore  Eq.  (33) 
should  constitute  a  good  benchmark  for  testing  discrete  scattering  kernels. 
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The  Diamond-Difference  SMART  Scattering  Matrix 

Our  two-dimensional  code  described  below  uses  diamond  differencing  in  both  space  and 
path  tength;  however,  only  the  numerical  treatment  of  the  path  length  variable  affects  the  definition 
of  the  scattering  matrix.  With  diamond  differencing  and  a  piecewise  constant  scattering  matrix,  it  is 
easy  to  show^*  that  the  solution  to  Eq.  (31)  is 

where  the  ${^1/3  are  the  edges,  AS{  the  size,  and  ^  the  SMART  scattering  matrix  for  the  i'th  path 
length  step. 

To  define  the  ^  we  force  Eq.  (36)  to  agree  with  the  exact  solution  [Eq.  (33)]  at  each  value 
of  It  is  shown  in  Ref.  15  that  the  required  scattering  matrix  is 

Si  -  ^  P  +  0(0-^., va)  G((>-*Si.i/3)-T* 

X  [G(0--Si^i/j)  G(0-.s,.i/j)-i-I]  .  (37) 

This  scattering  matrix  cancels  errors  due  to 
-angular  discretization 
-diamond  differencing  in  s 
-the  piecewise  constant  (in  s)  approximation. 

An  example  of  a  diamond  difference  SMART  scattering  matrix  is  given  in  Table  m  of  Ref.  15. 
Positivity 

We  demonstrate  in  Ref.  15  that  the  S'{  will  be  positive  provided  that  the  As^  are  sufficiently 
small,  and  the  basis  functions  [the  B'”(n)]  are  chosen  positive. 
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AltlMMigh  the  delta  function  corrected  scattering  matrix  is  SMART  if  Gauss  (luadrature 
points  are  used  (this  is  only  possible  in  one-dimensional  geometries)  it  b  not  positive.  (Cbmpare 
Table  V  with  Table  m  of  Ref.  IS.)  As  shown  in  Ref.  IS,  this  is  due  to  the  fact  that  the  basis 
functions  that  yield  the  delta  function  corrected  scattering  matrix  are  the  Lagrange  interpolation 
polynomials,  which  are  not  positive  functions.  With  positive  SMART  scattering  matrices  problems 
of  negative  fluxes  do  not  occur.  See  Fig.  1  of  Ref.  IS. 

SMART  scattering  matrices  have  produced  excellent  results  in  one-dimensional  single 
collision  Monte  Carlo^*  calculations  (see  Fig.  2  of  Ref.  IS)  in  two-dimensional  SR  calculations  (see 
Table  I  of  Ref.  17)  and  in  two-dimensional  S)f  calculations  (See  Figs.  3  and  4  of  Ref.  IS  and  Fig.  2 
of  Ref.  18). 

It  appears  that  SMART  scattering  theory  can  be  used  to  eliminate  or  reduce  several  of  the 
classic  problems  associated  with  electron  transport,  namely, 

-miniscule  mean-free-paths  that  imply 

a.  slowly  converging  source  iterations  in  So  and  SR  codes 

b.  many  collisions  per  history  in  Monte  Carlo  codes 

c.  relatively  large  ratios  of  cell-collided  to  cell-uncoUided  fluxes,  which  reduces 
the  accuracy  of  SR  calculations^ 

-large  angular  discretization  errors  in  Sjf  and  SR  codes 
-nonpositive  scattering  matrices  that 

a.  cannot  be  used  in  Monte  Carlo  codes 

b.  should  not  be  used  in  Sf^  codes  with  negative  flux  fixups 
-the  relatively  minor  problem  of  path  length  discretization  errors. 

SMART  scattering  matrices  have  performed  well  in  SR,  Si,;,  and  single  collision  Monte  Carlo 
computer  codes.  They  appear  to  be  an  attractive  alternative  to  Fokker-Planck  and  extended 
transport  correction  techniques  for  Sj.;  and  SR  calculations  and  to  the  condensed  history  approach 
for  Monte  Carlo  calculations. 
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Fokker-Hanck  owthods  give  a  good  representation  of  narrow  angle  scattering  but  are 
inap|»opriate  for  wide  angle  collisions. 

Extended  transport  corrected  cross  sections  can  simulate  both  narrow  and  wide  angle 
collisions,  and  are  in  fact  SMART  if  Gauss  quadrature  points  are  used.  Ftowever,  they  are 
nonpositive. 

With  the  condensed  history  approach,  direction  change  is  determined  by  sampling  the 
Goudsmit-Saunderson  distribution.  Since  this  distribution  is  the  solution  of  the  space-independent 
Spencer-Lewis  equation,  it  is  not  valid  near  material  interfaces  and  special  techniques  must  be  used.* 
As  speculated  above,  the  validity  of  a  SMART  scattering  matrix  should  not  depend  on  the  proximity 
of  a  boundary.  If  this  is  indeed  true,  then  the  SMART  scattering  matrices  should  prove  useful  in 
Monte  Carlo  codes  at  least  in  the  vicinity  of  interfaces.  To  date,  our  test  problems  have  involved 
only  media  with  vacuum  boundaries.  To  obtain  a  more  stringent  test  of  SMART  scattering  theory, 
calculations  with  several  closely  spaced  material  interfaces  will  soon  be  carried  out 

It  may  appear  that  the  method  for  generating  SMART  scattering  matrices  is  analogous  to  the 
flux- weighting  procedure  used  to  determine  effective  neutron  cross  sections.  Indeed,  both 
techniques  use  infinite  media  solutions  to  determine  appropriate  cross  sections.  However,  the 
SMART  scattering  matrix  is  by  no  means  an  angular  flux-weighted  average  of  the  true  scattering 
kernel.  In  fact,  the  total  cross  section  is  much  smaller  and  the  average  deflection  per  collision  much 
larger  than  those  for  the  true  scattering  kernel. 

SMART  scattering  theory  appears  to  have  an  important  advantage  over  flux-weighting 
techniques.  Flux-weighted  cross  sections  are  valid  only  for  one  particular  spectrum,  whereas  the 
SMART  scattering  matrix  gives  exact  solutions  for  any  arbitrary  angular  flux  (provided  it  is  spatially 
independent)  Thus,  flux-weighted  cross  sections  are  reliable  only  if  the  spectrum  in  the  region  to 
be  analyzed  resembles  that  used  to  generate  the  cross  sections.  SMART  scattering  matrices  should 
be  more  versatile.  Our  current  thinking  is  that  these  matrices  will  be  nearly  equivalent  to  the  true 
scattering  kernels  in  most  problem  geometries,  but  exactly  equivalent  only  in  infinite  media. 
Although  we  cannot  give  a  formal  proof,  our  reasons  for  this  assertion  are  as  follows: 
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-The  SMART  scattering  matrix  is  not  tied  to  a  particular  angular  flux  distribution. 

-There  is  considerable  numerical  evidence  that  the  SMART  scattering  matrices  work  well  in 
finite  media. 

-The  physical  model  corresponding  to  SMART  scattering  theory  (few  relatively  wide  angle 
collisions)  seems  to  be  a  reasonable  equivalent  to  the  true  situation  (many  small  angle  collisions). 

-Certainly  one  can  imagine  pathological  geometries  in  which  trajectories  with  a  large 
number  of  wnaii  angle  collisions  could  not  be  simulated  by  trajectories  with  a  lesser  number  of  large 
deflectioos.  Thus,  it  is  doubtful  in  our  opinion  that  an  exact  equivalence  can  be  demonstrated  for 
the  finite  medium  case. 

We  have  observed^*  that  the  SMART  and  conventional  scattering  matrices  are  in  close 
agreement  for  those  matrix  elements  that  represent  wide  angle  scattering.  However,  for  narrow 
deflections  the  elements  of  S  exceed  those  of  S.  The  augmented  values  of  these  elements  in  S 
compensate  for  the  large  number  of  very  narrow  deflections  that  are  too  small  for  direct  modeling 
on  the  Sn  quadrature  set 

The  conventional  scattering  matrix  fails  in  electron  transport  because  it  misses  the 
accumulative  effect  of  these  narrow  deflectioos. 

2.  SMART  First  Collision  Sources 

The  material  in  this  section  is  taken  from  Ref.  19.  The  use  of  analytic  first  collision  sources 
can  greatly  improve  the  accuracy  of  neutral  particle  S^  calculations.  For  electrons,  because  of  the 
strong  anisotropy  of  the  electron  scattering  kernel,  many  collisions  are  required  before  the  flux  from 
a  beam  source  is  smooth  enough  to  lend  itself  to  practical  numerical  treatment.  Thus  the  use  of  a 
conventional  analytic  first  collision  source  is  of  little  value. 

With  SMART  scattering  theory^*  as  described  above,  the  true  electron  scattering  law  is 
replaced  by  an  equivalent  one  that  involves  many  times  fewer  collisions  but  larger  angular 
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d^toctioas  per  collision.  With  such  a  scattering  law,  the  electron  transport  emulates  that  of  neutral 
particles.  Hence,  application  of  an  analytic  first  collision  source  should  be  effective.  Ifowever, 
^lART  scattering  cross  sections  are  defined  in  terms  of  a  scattering  matrix  and  thus  only  apply  for 
scatter  from  one  quadrature  direction  to  another. 

We  have  extended  the  theory  to  include  scatter  from  (but  not  to)  an  arbitrary  beam 
direction.  This  entails  the  generation  of  one  additkinal  column  of  the  SMART  scattering  matrix. 
The  m'th  element  ^  of  the  vector  is  a  SMART  cross  section  for  transfer  from  the  beam  to 

A 

direction  flP. 

The  ^  are  determined  by  a  procedure  similar  to  that  used  to  determine  the  SMART 
scattering  matrix.  We  require  that  the  Sjf  method  with  an  analytic  first  collision  source  produce  the 
exact  solution  to  an  appropriate  benchmark  problem. 

Again  we  choose  Eq.  (31)  to  define  our  benchmark  problem  and  begin  with  the  discrete 
ordinates  solution.  We  assume  that  the  SMART  scattering  matrix  S  has  already  been  determined. 
Then,  for  the  fast  path  length  step  it  is  easy  to  show^*  that  the  S{f  expression  for  the  multiply 
scattered  flux  is: 


7'(As)  -  jl  -  ^sj  (38) 

where 

^  si 

m*l 

Wq,  «  weight  for  direction,  CP* 

m  vector  whose  m'th  component  is  the  multiply  scattered  flux  at  (7”  and  s  «  As.  Let 


the  exact  (continuous  angle)  expression  for  the  multiply  scattered  flux  be  denoted  byT(^)-  Then, 
the  exact  and  the  S{  solution  agree  provided  that  we  pot 
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lb 


1-  exp(-ab 


aodl‘(As)  must  satisfy  the  normalizatioa^* 


(40) 


S 

m^l 


"m  8S,(As)  -  1  -  exp(-Ub  As) 


(41) 


These  results  require  careful  interpretation.  Equations  (40)  and  (41)  are  insufficient  to 
determine  both  ami  because  the  maximum  component(s)  ofl‘(As)  is  (are)  not  defined.  Here 
T^(As)  represents  the  multiscattered  rather  than  the  total  flux.  Ordinarily,  these  two  fluxes  differ 
only  by  a  delta  function  of  negligible  strength  [exp(-oAs),  where  o  is  the  true  cross  section]  centered 
at  the  beam  direction.  However,  as  explained  above,  it  is  necessary  that  oi,«  o.  This  amounts  to 
broadening  the  definition  of  the  uncollided  beam  to  include  some  narrowly  scattered  electrons, 
thereby  reducing  the  maximum  element(s)  of  l<‘(As). 

The  amount  of  broadeiUng  is  determined  by  which  remains  a  free  parameter.  Once 
is  chosen,  the  maximum  element(s)  of7^(As)  can  be  determined  from  normalization  Eq.  (41)  and 
from  Eq.  (40).  This  choice  of  guarantees  that  the  ^  method  will  reproduce  the  exact 
benchmark  result,  that  is,  the  normalized  1*  (As). 

Unlike  C),,  the  attenuation  facton  for  travel  along  any  of  the  quadrature  directions  are 
uniquely  determined  by  SMART  scattering  theory,**  and  they  all  have  comparable  values.  Our  Sjif 
code  selects  one  of  these  (corresponding  to  m  a  k,  for  example)  for  O),.  Because  of  this  choice  and 
the  corresponding  broadened  meaning  of  the  uncollided  flux,  the  beam  is  not  modeled  by  a 
monodirectional  ray;  imtead  it  is  modeled  by  a  cluster  of  rays  with  an  angular  spread  commensurate 
with 

Figure  la  of  Ref.  19  shows  the  result  of  the  application  of  the  Siif  method  with  the  SMART 
first  collision  source  to  a  two-dimensional  problem.  A  beam  of  200-KeV  electrons  is  assumed 
incident  on  a  0.01  x  0.02  g/cm*  aluminum  slab.  The  beam  obliquity  (in  the  x-y  plane)  is  45  deg. 
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and  die  point  of  ii^idence  is  located  at  x  >  0,  y  •  0.01  (as  shown  on  the  inset).  The  rectangular 
slab  is  divided  into  50  square  zoms  (0.002  x  0.002),  and  the  histogram  displays  the  energy  (in  kilo- 
electron-volts)  deposited  in  each  zone  as  obtained  using  an  calculation.  Hgure  lb  of  Ref.  19 
shows  the  result  of  the  application  of  the  ACCEPT  (Ref.  20)  Monte  Carto  code  to  the  same 
problem.  The  two  methods  agree  welL  The  S4  calculation  consumed  ^'th  of  the  computer  time 
required  by  the  Nfonte  Carlo  code  (100,000  case  histories).  It  is  anticipated  that  better  agreement 
would  be  achieved  with  a  higher  order  calculation  such  as  S,,  which  would  run  approximately  three 
to  four  times  slower  than  the  S4. 

B  Extension  of  the  Streaming  Ray  (SR)  Method  To  Two-Dimensional  Homogeneous  Media 

This  section  is  condensed  from  Ref.  17.  The  method  of  streaming  rays  (SRs)  has  been 
shown  previously^  to  be  an  elective  algorithm  for  one-dimensional  electron  transport  studies.  The 
Spencer-Lewis  equation*"'^  is  solved  for  the  electron  distribution  in  direction  fi,  position  x,  and  path 
length  s,  and  the  continuous  slowing  down  approximation  is  used  to  relate  energy  loss  to  distance 
traveled,  so  that  path-length-dependent  cross  sections  can  be  defined.  The  SR  algorithm  for 
electron  transport  calculations  is  advantageous  because  it  facilitates  modeling  of  the  continuous 
energy  loss. 

We  have  extended  the  SR  method  to  two  spatial  dimensions  (three  phase-space  dimensions, 
X,  y,  and  s). 

In  contrast  to  its  one-dimensional  predecessor,*  the  SR2D  code  accommodates  nonuniform 
cell  dimensions  in  x  and  y  and  allows  for  arbitrary  discrete  ordinates  quadrature  sets  (S,,  S4,  Sg,  S,, 
S^,  or  Sjg).  Families  of  streaming  rays  originate  in  the  x  >  0  plane  with  a  uniform  spacing  and 
overlay  the  three-dimensional  Eulerian  grid  in  x,y4  phase  space.  Each  ray  is  defined  by  its 
direction  m  and  the  coordinates  of  its  origination  point  in  the  y,s  plane.  With  this  arrangement  and 
with  As,  the  path  tength  step  equal  to  an  integer  multiple  of  DS,  the  spacing  between  streaming  ray 
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originatioii  points,  each  slab  of  thickness  ^  has  an  identical  pattern  of  interlacing  streaming  rays 
and  path  lengths  need  to  be  computed  only  once. 

The  SR2D  code  uses  SMART  cross  sections^*  that  significantly  reduce  angular  dttcretization 
erron.  This  allows  the  highly  anisotropic  electron  scattering  to  be  modeled  with  relatively  few 
discrete  directions. 

The  code  was  used  to  calculate  the  energy  deposition  profile  for  an  isotropic  point  source  at 
the  periphery  of  a  two-dimensional  aluminum  medium  with  dimensions  0.01  g/cm>  thick  x  0.02 
g/cm*  wide.  The  computational  grid  was  5  x  10  uniform  cells,  respectively.  The  path-length 
increment  was  0.002  g/cm*  with  25  path-length  increments  chosen.  An  S,  qxiadrature  set  was  used. 
The  isotropic  point  source  was  normalized  to  one  incident  particle  with  an  energy  of  200  KeV. 

To  validate  the  SR2D  resxilts,  the  test  problem  was  also  solved  using  the  electron/photon 
Monte  Carlo  code  TIGER  (Ref.  9).  The  total  energy  deposited  in  the  medium  and  peak  cell 
energies  was  selected  to  facilitate  the  comparison.  Results  for  SR2D  and  TIGER  are  provided  in 
Table  I  of  Ref.  17.  For  this  problem,  the  calculated  values  of  the  total  energy  deposited  in  the 
aluminum  were  within  1%,  but  peak  cell  energies  varied  by  4%.  The  largest  relative  error  was 
<30%,  and  this  occurred  where  numerical  values  were  small,  well  away  from  the  area  of  peak 
energy  deposition.  The  TIGER  results  came  from  an  evaluation  of  50,000  case  histories  and  were 
within  ±9%.  Other  problems  have  been  compared,  particularly  a  monodirectional  point  source,  with 
equally  good  results. 

C.  The  Two-Dimensional  Multiregion  /Diamond  Difference  Algorithm 

Although  the  SR2D  code  gives  good  results,  the  algorithm  would  be  significantly  more 
complicated  for  multiregion  problems.  We  therefore  decided  to  use  initially  the  somewhat  less 
complicated  method  for  our  multiregion  code. 
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I.  Theory 

Extei^ing  from  homogeneous  to  multiregions  introduces  a  complication.  In  the  Spencer- 
Lewis  equation  s  is  used  to  specify  the  electron's  energy.  For  homogeneous  media,  the  CSDA 
implies  a  one-to-one  correspondence  between  E  and  s.  However,  with  multiregion  problems  this 
one-to-one  correspondence  is  lost  (since  the  stopping  power  is  spatially  dependent)  and  s  no  longer 
determines  E. 

Therefore,  we  redefine  s  to  mean  the  path  length  that  would  be  required  to  reach  energy  E 
in  an  infinite  homogeneous  medium  of  the  composition  occurring  at  the  electron's  location.  That  is, 
we  let 


s  ■ 


dE' 

f  (E') 


(42) 


where 


ds 


(1) 


is  the  stopping  power  in  region  r.  Thus  s  is  not  the  true  path  length  for  an 


electron  that  has  traveled  in  more  than  one  region. 

With  this  definition,  E  is  determined  by  s  and  the  electrons  location.  However,  the  electron 
flux  at  a  particular  value  of  s  is  no  longer  continuous  across  mterial  interfaces.  Thus,  the 
equations  cannot  be  formulated  in  the  conventional  maimer.  However,  the  flux  integrated  over  a 
path  length  step  will  be  continuous  provided  that  the  path  length  steps  represent  the  same  energy 
intervals  in  each  region.  By  substituting  integrated  quantities  for  cell  center  fluxes  and  some  of  the 
edge  fluxes  we  are  able  to  use  the  standard  algorithm. 

We  consider  first  the  standard  discrete  ordinates  from  of  Eq.  (1)  in  x-y-s  geometry,  for 


mesh  ceil  Axj  Ay^  in  region  n 


(43) 


In  this  equation  whole  indices  i,j.k  are  used  for  quantities  that  have  been  averaged  over  Asf,  Axj 

and  Ay](  respectively,  while  half  indices  represent  cell  edge  values.  (The  need  for  the  superscript  r 

on  is  explained  below).  To  better  model  the  extreme  anisotropy  of  the  scattering  kernel  the 

are  determined  using  SMART  scattering  theory. 

The  Su  algorithm  cannot  be  applied  directly  to  £q.  (43)  because  the  edge  fluxes  ,  and 

ij±|k 

,  are  not  continuous  at  material  interfaces.  However,  these  same  edge  fluxes  when  .multiplied 
ijktj 

by  Asf  are  continuous  provided  that  the  Asf  define  the  same  energy  Interval  in  each  region,  that  is, 
provided  that 


where  the  energies  £.  ^  defining  the  edges  of  the  i'th  path  length  step  are  identical  in  each 
**  2 

material  region.  Multiplying  Eq.  (43)  by  As'  we  obtain 
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Ay^ 


Uk+i 


SLn'^  +  QS 


m'-l 


(45) 


where  the  bar  (-)  is  used  to  indicate  quantities  that  have  been  multiplied  by  As'.  All  edge  fluxes  in 
£q.  (45)  are  continuous,  and  except  for  the  numerical  treatment  of  the  scattering  term,  the  equation 
is  exact.  However  there  are  four  unknowns,  the  cell  center  flux  and  the  three  cell  exit  fluxes.  In 
order  to  obtain  a  solvable  system  of  equations  we  use  the  diamond  difference  approximation,^^ 
which  here  takes  the  form 


(46) 


Except  for  the  appearance  of  Asf  in  Eq.  (46)  and  its  absence  from  Eq.  (45),  these  last  two 
equations  are  identical  in  form  to  the  conventional  /Diamond  Difference  equations.^**^*'  Making 
a  minor  adjustment  for  Asf  we  solve  these  equations  using  the  usual  Si,( /Diamond  Difference 
algorithm. 


2.  Results 

We  have  performed  many  two-dimensional  multiregion  problems.  The  results  are 
qiialitatively  correct  and  the  calculations  conserve  particles  to  at  least  four  significant  figures.  This 
is  encouraging,  however  we  have  not  yet  compared  the  multiregion  results  to  accurate  benchmark 


solutions. 


-  23  - 


The  code  has  been  verified  with  several  single  region  problems.  IHgures  3  and  4  and  Table 
vn  of  Ref.  IS  show  a  comparison  of  and  Monte  Carlo  calculations  of  the  energy  deposition 
profile  due  to  a  200  KeV  isotropic  point  source  of  electrons  incident  on  a  two-dimensional 
aluminum  slab.  Similar  comparisons  for  both  isotropic  and  monodirectional  beam  sources  are  shown 
in  Figs.  1  and  2  of  Ref.  18.  The  two  methods  are  in  very  good  agreement 


3.  Description  of  the  SN2D  Code 

a.  Input 

The  user  may  select 

-Number  of  regions  (1  to  15) 

-Number  of  elements  per  region  (1  to  10) 

-Number  of  energy/path  length  mesh  ceils 

-Number  of  x-mesh  cells 

-Nmnber  of  y-mesh  cells 

-Sjif  number  (2,  4,  6,  8,  12  or  16) 

-Input  in  centimeters  or  micrometers 

-Maximum  number  of  source  iterations  per  path  length  step 

-Reduced  or  full  printout  option 

-Type  of  source 

1.  Spatially  distributed  isotropic 

2.  Spatially  distributed  monodirectional  along  one  of  the  quadrature 

directions 

3.  Point  source  at  user  specified  location  in  user  specified  direction 
-Energy  of  source  particles 

-Negative  flux  fix-up  option 
1.  No  fix-up 
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2.  Negative  flux  fix-up 

3.  Negative  flux  fix-up  and  reduce  As  until  scattering  matrix  is  positive 

4.  Negative  flux  and  negative  source  fix-up 

5.  Start  with  2  and  convert  to  1  if  the  scattering  matrix  goea  negative  at 
wme  path  length  step  (we  recommend  option  3) 

-Beam  source  representation 

1.  Monodirection 

2.  Ouster  of  four  nearly  monodirectional  beams.  (This  better  represents  the 
slight  spreading  of  a  SMART  "uncollided"  beam  as  explained  in  Section 
n.A2 

-Cross  section  option 

1.  Screened  Rutherford 

2.  Riley  (not  fully  debugged) 

-Maximum  order  of  Legendre  cross  section  expansion  for  determining  the  Goudsmit- 
Saunderson  matrix.  (The  suggested  value  is  200) 

-The  energy  mesh  intervals  (arbitrary  spacing) 

-The  x-mesh  intervals  (arbitrary  spacing) 

-The  y-mesh  intervals  (arbitrary  spacing) 
b.  Output 

The  output  from  the  code  includes 
-Angular  fluxes  for  each  cell 

-The  spatially  integrated  angular  flux  in  each  energy  step 
-The  total  flux  for  each  mesh  cell 

-The  total  leakage  currents  (electrons  per  second)  from  each  surface 
-The  energy  deposited  in  each  spacial  mesh  cell. 
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D.  The  Super  Computer  (SCSn)  Method 

The  material  in  this  section  is  taken  from  Ref.  21.  When  devetoped  in  the  'SOs  and  '60s,  the 
discrete  ordinates  method’^  was  tailored  for  optimum  efficiency  on  the  computers  of  the  day. 
Computation  speeds  and  the  available  memory  of  those  fint  computers  were  extremely  poor  by 
today's  standards.  Nevertheless,  the  algorithm  was  and  continues  to  be  the  method  of  choice  for 
many  transport  problems.  However,  Si^  calculations  have  been  plagued  by  two  serious  problems, 
numerical  diffusion  and  ray  effects;  the  first  arises  from  spatial  discretization,  and  can  be  quite 
severe^  near  flux  discontinuities;  the  second  results  from  angular  discretization  and  is  very 
pronounced  for  low  scattering  media.  Although  the  SR  method^*^^  eliminates  the  former  problem  it 
does  not  eliminate  the  latter.  We  have  developed  a  new  discrete  ordinates  method  specifically 
designed  for  todays  large  memory,  vector  and/or  parallel  processing  computers  that  sharply  reduces 
numerical  diffusion  and  eliminates  observable  ray  effects. 

The  method  has  been  tested  with  excellent  results  on  one>dimension  homogeneous  slab 
problems.  For  these  problems  the  Spencer  Lewis  equation  [£q.  (1)]  for  discrete  direction 
becomes 


^x,;*“,s)  -  J  S^'  +  Q(x,m",s)  ,  (47) 

m'-l 

where 

Discretizing  path  length  and  space  yields  the  foUowing  equation  that  is  solvable  by  the 
standard  source  iteration  method, 


where 


(49) 


-  26  - 


w7«'(p)+'^«  .  (50) 

m'-l 

and  ^  have  components  and  Q(Xj,m'^^),  j  and  i  are  the  x  and  s  discretization 

indices,  and  Lm  now  the  matrix  representation  of  the  original  operator.  It  is  evident  from  Eq.  (3) 
that  each  matrix  element  of  (Lni)~^  represents  the  uncollided  flux  in  a  particular  cell  due  to  a 
monodirectional  unit  source  in  another.  Therefore  is  a  very  large  and  sparse  matrix.  The 

conventional  algorithm  enables  Eq.  (49)  to  be  implemented  without  explicitly  determining  the 
matrix  elements  of  (L^,)'^,  or  explicitly  executing  the  implied  multiplications  by  zero.  This  is  an 
enormous  advantage  for  low  memory  scalar  computers,  but  it  is  also  the  major  source  of  the 
numerical  errors.  Although  the  matrix  elements  corresponding  to  adjacent  cells  are  determined 
accurately  (implicitly),  those  corresponding  to  distant  cells  are  not,  because  they  are  obtained 
(implicitly)  through  the  repeated  application  of  a  spatial  differencing  approximations.  Some 
elements  which  should  be  zero  turn  out  positive  (or  even  negative  if  no  fix-up  is  used)  resulting  in 
numerical  diffusion.  Matrix  elements  for  cells  lying  along  a  discrete  direction  tend  to  be  over¬ 
estimated  causing  ray  effects.  Since  scattering  tends  to  reduce  the  magnitude  of  the  matrix  elements 
for  the  distant  cells  it  also  reduces  the  severity  of  these  two  errors. 

The  fundamental  ideas  of  our  method  are  simple: 

-Calculate  the  elements  of  (L^,)*^  explicitly  and  accurately. 

-Implement  Eq.  (49)  in  matrix  form. 

Item  2  above  allows  us  to  take  advantage  of  vector  and  or  parallel  processing  and  thereby 
compensate  for  the  many  multiplications  by  zero. 

Since  (L„)~^  represents  an  operator  for  uncoUided  transport,  its  elements  can  be  calculated 
to  any  degree  of  accuracy.  However,  this  would  require  the  evaluation  of  tedious  integrals,  and 
instead  we  have  adopted  the  following  simple  scheme. 

The  first  step  is  to  calculate  the  uncollided  flux  contribution  from  a  monodirectional  unit 
source  in  cell  i',  j'  to  the  target  path  length  cell  i  by: 


1  1.0  - 

ffj  '  6s^•(a^y*  ’  * "  *  ’ 

4  -  (51) 

AS{  (7]'  O’]  ,  -  -  -  . 

Although  our  code  is  more  general,  for  simplicity  in  £q.  (51)  we  have  assumed  that  the 
mesh  cells  are  square.  Both  expressions  are  derived  by  performing  an  uncollided  electron  balance 
over  the  i'th  path  length  cell.  The  uncollided  flux  is  then  distributed  to  the  individual  spatial  cells 
in  proportion  to  the  area  subtended  on  the  target  cell  by  an  angular  cone,  originating  from  the 
source  cell 

For  this  angular  cone  we  use  one  of  two  forms,  depending  on  the  relative  distance  of  the 
target  path  length  from  the  path  length  of  the  source  celL 

Figure  1  illustrates  the  constniction  of  the  angular  cone  for  the  discrete  direction  m,  used 
for  target  path  lengths  when  i  >  i'  1.  The  angular  cone  is  generated  by  extending  the  region  of 
influence  of  the  discrete  direction  m  to  include  half  of  the  angular  distance  to  each  of  the  m-fl  and 
m-1  directions.  The  area  created  by  the  intersection  of  this  angular  cone  and  the  i  ±  |  bondaries 
of  the  target  path  length  cell  defines  a  region  over  which  is  to  be  distributed.  The  individual 
elements  of  are  then  calculated  based  on  the  ratio  of  area  contained  within  a  particular  mesh 

cell  to  this  total  area. 


where 

Area  y  «  area  of  intersection  within  cell  ij 
and 

Area  ]  -  total  area  of  intersection  within  path  length  cell  i 
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If  a  ceU  within  the  target  path  length  is  not  intercepted  by  the  angular  cone  the  matrix  element  for 
that  target  cell  is  lero.  In  fact  the  majority  of  these  elements  are  zero,  as  is  clearly  evident  from 
Fig.  1. 

Figure  2  illustrates  the- construction  of  the  angular  cone  for  target  path  length  cells  for 
which  i  ^  i'  4-  1.  Whereas  tlw  angular  cone  originated  from  the  center  of  the  source  cell  for  the 
target  path  length  steps  above,  it  is  now  constructed  so  that  the  m  ±  ^  boundaries  of  the  cone, 
originate  from  the  centers  of  the  spatial  edges  of  the  source  cell  Once  the  area  is  defined  £q.  (52) 
is  again  used  for  calculating  the  individual  elements. 

The  distinction  between  these  two  cases  is  made  based  on  physical  arguments.  If  the 
previous  cell  centered  scheme  were  used  throughout  the  discrete  ordinates  grid,  the  elements  of 
[Lm]-i  for  adjacent  cells  would  be  underestimated.  In  fact,  a  simple  examination  of  Fig.  1  shows 
that  for  this  scheme,  regardless  of  the  discrete  direction,  the  elements  of  [Lg,]~^  for  the  cells 
immediately  above  and  below  the  source  cell,  i'j',  would  be  zero. 

On  the  other  hand,  if  the  above  cell  edge  center  scheme  (Fig.  2)  were  applied  to  the  entire 
grid,  the  extended  regions  of  adjacent  discrete  directions  would  overlap  for  distant  cells 
overcompensating  for  the  use  of  discrete  ordinates. 

This  semi-analytic  geometric  scheme  described  above  for  computing  the  elements  is 

the  fundamental  mechanism  whereby  both  ray  effects  and  numerical  diffusion  are  removed.  The 
use  of  discrete  ordinates  is  compensated  for  by  extending  the  region  of  influence  of  each  discrete 
direction  to  include  a  portion  of  the  previously  'unseen'  space  between  adjacent  directions.  In  this 
manner  the  streaming  particles  are  allowed  access  to  the  entire  phase  space.  The  numerical 
diffusion  has  been  reduced  by  explicitly  setting  to  zero  those  elements  of  [Lnr^  which  should  be 
zero. 

The  use  of  cones  (instead  of  discrete  directions)  is  not  a  new  idea,*^  however  here  they  are 
used  to  determine  all  of  the  matrix  elements  directly,  and  not  just  those  corresponding  to  cells 
inunediately  adjacent  to  the  source  celL  It  is  the  use  of  cones  for  uncollided  transport  to  distant 
cells  that  eliminates  ray  effects.  Unlike  other  ray  effect  mitigation  techniques  this  new  method  is 
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quite  effident,  being  2.6  times  slower  than  conventional  Sn  on  scalar  machines,  but  only  1.2  times 
^wer  on  vector  machines. 

Because  of  the  many  zeros  in  the  matrix  the  SCS^  method  does  not  make  efficient 

use  of  computer  menmry  nor  does  it  make  efficient  use  of  the  computational  resources  of  a  scalar 
(due  to  the  many  unnecessary  multiplicatioia  by  zero).  However  with  the  large  memory 
nwtfhiiMa  now  available  the  former  inefficiency  is  often  irrelevant  and  the  fact  that  the  source 
iteration  ‘a  cast  in  matrix  form  enables  extensive  use  of  vector  and/or  parallel  processing  that  can 
compensate  for  the  latter  inefficiency. 

The  SCSn  method  was  tested  using  a  200  KeV  isotropic  point  source  incident  on  aluminum 
slabs.  With  coarse  mesh  cells  (37  mfp  squares)  the  new  and  conventional  methods  were  in  close 
agreement  (see  Fig.  3).  However,  with  fine  mesh  cells  (3.7  mfp  squares)  the  new  method  was  not 
subject  to  the  ray  effects  evident  with  the  conventional  solution  (see  Fig.  4). 

To  demonstrate  the  lack  of  numerical  diffusion  obtained  with  the  SCSf(  method  we  have 
calculated  the  electron  flux  due  to  a  monodirecdonal  source  ip  •  .30)  incident  on  the  s  •  0  surface 
of  a  vacuum.  The  results  are  shown  in  Table  I  along  with  those  obtained  from  the  conventional 
method  with  and  without  a  negative  flux  fix-up,  and  in  Table  n  along  with  those  obtained  using 
the  SR  method  with  and  without  a  ray  effect  mitigation  routine. 

The  SCS;f  results  are  clearly  superior  to  either  of  the  conventional  solutions  in 
representing  the  discontinuity  along  the  edge  of  the  beam.  If  the  source  particles  are  interpreted  as 
a  nmnodirectional  beam,  then  the  SR  results  with  no  ray  effect  mitigation  represent  the  true  solution 
and  all  particles  should  remain  within  the  dotted  lines  shown  in  Table  n.  On  the  other  hand,  if  the 
source  particles  are  assumed  to  be  distributed  between  p  -  25  and  p  •  .35,  the  edges  of  the  cone 
corresponding  to  direction  n  •  .30,  then  the  flux  should  be  contained  between  the  solid  lines  shown 
in  Tables  I  and  n,  and  the  results  are  best  calculated  by  the  SCS^  method.  For  this  interpretation 
of  the  source,  the  SR  method  does  better  with  the  ray  effect  mitigation  routine  turned  on,  however, 
it  is  evident  that  this  routine  is  only  partially  successful,  as  the  flux  still  peaks  along  the  source 


direction. 
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HL  Numerical  Diffusion  Theory 

The  mathraiatical/numerical  descriptioa  of  how  electrons  interact  with  electronic  devices  can 
be  classified  into  three  fundamental  calculational  categories.  The  probabilistic  approach  embodied  in 
the  Monte  Carlo  method  is  simply  a  numerical  simulation  of  the  collision  interactions  and 
subsequent  streaming  of  electrons  and  requires  considerable  computational  resources  to  use 
e^bcdvely.  The  deterministic  transport  methods  as  described  in  Sections  the  previous  sections  offer 
a  viable  alternative  to  the  Monte  Carlo  method  but  may,  in  some  cases,  provide  more  detailed 
information  than  is  necessary  for  an  adequate  description.  In  particular,  at  the  low  end  of  the 
electron  energy  spectrum  when  the  electrons  have  experienced  many  collisions,  the  complete  angular 
flux  variation  in  angle  is  not  always  necessary  since  the  flux  is  nearly  isotropic.  For  this  reason,  an 
effort  has  been  initiated  to  develop  a  diffusion  theory  description  in  parallel  with  the  transport 
theory  description.  The  diffusion  equation  is  a  macroscopic  deterministic  description  of  electron 
motion  and  involves  far  less  computational  effort  than  either  Monte  Carlo  or  transport  theory 
methods. 

The  main  goal  of  the  diffusion  theory  component  is  to  generate  a  reliable  three-dimensional 
algorithm  that  will  be  coupled  to  a  2-D  or  3-D  algorithm  at  an  appropriately  low  energy.  In 
this  way,  we  hope  to  treat  both  the  high  and  low  energy  regions  with  the  greatest  efficiency. 

In  the  following,  the  progression  toward  the  completion  of  the  3-D  algorithm  is  documented. 
We  begin  with  a  1-D  development  which  establishes  an  adequate  foundation  for  the 
multidimensional  algorithms  to  follow.  From  this  initial  study,  a  proper  undentanding  of  how  the 
multidimensional  algorithm  is  obtained.  Many  of  the  numerical  features  of  the  1-D  treatment  have 
been  incorporated  into  the  multidimensional  algorithms  making  the  1-D  study  invaluable.  An 
alternative  approach  would  have  been  to  use  an  "off-the-shelf  diffusion  code  with  enough 
flexibility  to  accommodate  electron  motion  in  condensed  matter;  however,  the  benefit  of  the  very 
worthwhile  learning  experience  would  have  been  lost 
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A.  1-D  DifTuskm  Theory 

The  generic  diffusion  equation  describing  electron  transport  as  modeled  by  the  Spencer* 
Lewis  equation  (1)  can  be  written  in  one-dimension  as 

^  (X4)  -  ^  P(X4)  ^  -  q(x4)  ^X4)  +  S(x^)  (53a) 

where 

a  electron  angularly  integrated  density  (simply  called  density)  for  electrons  having 
transversed  a  total  path  length  s  to  arrive  at  position  x 

p  a  diffusion  coefficient 

q  a  absorption  cross  section 

S  a  external  source. 

The  density  is  subject  to  the  following  initial  and  general  mixed  boundary  conditions 

^x,0)  a  f(s)  (53b) 

^  -  o(x)  at  X  a  a,  b  .  (53c) 

at  the  boundaries  a  and  b  where  7^  are  known  functions.  Equations  (53)  can  be  obtained  from 
Eq.  (1)  in  many  ways.  The  simplest  method  is  to  integrate  Eq.  (1)  over  n  (in  one  dimension)  and 
define  the  density  as 


^x^) 


I  ^dM  0(x,/*4) 


and  current  as 


J(X4)  ■ 


dp  p  ^xpi.s) 
-I 
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and  then  invoke  Pick's  law 

J(x4)  “  -  ^  P(*^)  ^  ^**»)  • 

The  absorption  term  is  included  for  generality  and  the  boundary  conditions  result  from  partial 
crinent  considerations.  Implicit  in  the  derivation  of  Eqs.  (53)  is  the  assumption  of  a  linear  angular 
dependence  of  the  angular  flux  on  the  electron  direction  ft  in  one  dimension.  This  is  true  only  after 
many  collisions  and  a  distance  of  at  least  three  mean  free  paths  from  discontinuities  such  as  material 
boundaries  and  spatially  localized  sources.  For  tlmse  reasons,  Eqs.  (53)  are  anticipated  to  be 
applicable  when  the  electrons  have  traversed  a  relatively  large  path  length  at  (low  energy)  and  in 
large  homogeneous  regions.  Nevertheless,  a  reasonable  approximation  can  be  obtained  from 
diffusion  theory  at  material  boundaries  of  scattering  materials. 

1.  Poiat  Spadil  Finite  Difference  Scheme*^ 

To  derive  a  numerical  algorithm  to  solve  Eqs.  (S3),  we  first  partition  the  x  variable  into 
discrete  intervals  as  shown  in  Fig.  5.  The  interval  midpoint  is  designated  as  i  and  the  cell  right  and 
left  boundaries  by  i-  ^  and  i^t-  ^  respectively.  Next  Eq.  (S3a)  is  integrated  over  cell  i  to  give 


dx  ^x,s)  - 


g<«>- 


dx  q(x,s)  0(x4)  + 


dx  S(x,s) 


Then  integrating  the  second  term  exactly  and  using  the  approximations 
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I  *dx  h(x4)  ss  ^  ^  ^(s) 


dx 

» (^(s)  -  ^-i(s))/AXi 

Mxji) 

ax 

-  -  0|(s))/Ax,+i  , 

*■  1 

we  arrive  at  the  spatially  differenced  equation 

5*<») 

-  %(«)  +  Ci(s)  [^^.iCs)  -  ^(s)]  -  r,(s)  ^(s)  +  Q|(s) 

(54a) 

where 

■ 

2p,(s) 

Axi(Axi+AXi+i) 

(54b) 

Ci(s)  ■ 

2  Pi*i(s) 

Ax,+i(Axi+Axj+i) 

(54c) 

ri(s)  - 

(qi(s)  AXi  +  qi(s)  Axi+i)/(Axi+Axt+j) 

(54d) 

Qi  - 

(St(s)  Axi  +  S,+i(s)  AXj+i)/(AXi+AXi+i)  . 

(54e) 

At  this  point  there  are  several  ways  the  path  variable  s  could  be  handled  and  in  the  following 


sections  we  discuss  two  of  these. 
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2.  Continuous  Analytical  Continuation 

A  theoretkaily  simple  and  computationally  efficient  method  for  effecting  the  s  integration 
of  Eqs.  (54)  is  found  by  using  Taylor  series.**  This  method  has  been  rediscovered  recently  and  it  is 
quite  surprising  that  it  has  not  found  more  widespread  use. 

To  begin  with,  we  assume  that  P{(s),  q|(s)  and  Si(s)  are  given  functions  of  s  in  a 
homogeneous  cell  i  and  these  functions  can  be  expanded  in  Taylor  series  in  s,.!  <  s  <  s,  to  give 


w  r 

a,(5).  ^  «5a) 

k-0 


Ci(s) 

OO 

-I 

(55b) 

k-0 

rj(s) 

OO 

- 1 

(55c) 

k-0 

Qi(s) 

OO 

- 1 

(55d) 

k-0 

According  to  Eqs.  (S4b-e)  the  expansion  coefficients  Cy^,  and  are  related  to  the 
expansion  coefficients  of  pj,  q{  and  Sj  around  an  arbitrary  point  s^.^.  When  ^(s)  is  similarly 
expanded  in  Sf.i  <  s  <  s, 

OO  ,f 

^(s)-  ^  .  (56) 

k-0 

a  recurrence  relation  results  for  when  Eq.  (56)  is  introduced  into  Eq.  (S4a) 
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“  S  (iS)!  ■ 

j-0 


k! 


j-0 


(k-j)! 


Cf*-i 


i+lj 


-  ^j]  -  ^  Itk-j  ^  ■*■  ^Jt-1 
j-0 


with 


^r.o  -  ^(vi) 


(57a) 


(57b) 


Thus,  we  have  an  explicit  two-dimensional  recurrence  relation  for  ^  which  when  introduced  into 
Eq.  (56)  provides  the  solution  ^(s).  The  advantage  of  this  formulation  is  that  round  off  error  can 
be  controlled  by  the  choice  of  the  interval  boundaries  s,.  These  points  are  chosen  such  that  the 
series  in  Eq.  (56)  is  convergent  and  requires  less  than  about  ten  terms  for  convergence  to  a  desired 
accuracy.  By  using  the  previously  calculated  value  ^(s^i)  as  the  initial  condition  in  the  interval 
Vi  ^  s  5  s,,  we  are  continually  analytically  continuing  the  Taylor  series  representation  to  the 
subsequent  intervals— hence  the  name  of  the  method.  This  approach  has  been  shown  to  be  very 
accurate  in  other  areas  of  application  such  as  neutron  slowing  down  theory,  extended  gas  kinetics 
and  radioactive  decay.  A  second  major  advantage  of  the  method  (not  used  here)  is  in  application  to 
nonlinear  ordinary  differential  equations. 

The  method  can  equally  be  applied  to  the  multidimensional  diffusion  equation  and  its 
numerical  implementation  will  be  disctissed  for  the  3-D  application  in  a  future  report. 

3.  Path  Length  Discretization 

Following  along  more  conventional  lines,  the  path  length  derivative  can  be  descretized  as 


As 


-  ^")/As 


follows: 
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In  addition,  if  ^(s)  is  composed  of  fractions  of  and  in  the  form 

^(s)  -  o  4^*^  +  (1-a)  0  <  o  <  I, 

Eq.  (S4a)  becomes 

[i  -  a  6r]  A"**  -  ♦  (l-o)  5."  »!"  ♦  s]  * 


+  (1-a)  cj*  +  ^  (58a) 

where 

if  •  (af  +  af"^)/2  (58b) 

cf  n  (cf  +  cf"^)/2  (58c) 

ff  •  (rf  +  rf*'^)/2  (58d) 

bf  •  -ii-Ci-fi  (58e) 

^  >  (Qf  +  Qf"^)/2  .  (580 


In  matrix  form,  Eq.  (S8a)  can  be  written  as 

.  a  A"  +  (1-a)  A“  7"  +  Q" 

A«  —  — 


^  -  ^(Sn) 


(59) 
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where  A”  is  the  classic  tridiagonal  matrix  which  will  be  exploited  below. 

By  varying  the  value  of  a,  one  obtains  the  following  well-known  discretization  schemes 
a «  1  ,  implicit 
a «  0  ,  explicit 
a «  I  ,  Crank-Nicolson. 


4.  Numerical  Implementation  and  Test  Case  Results 
a.  Iteration 

Iteration  is  commonly  used  to  solve  Eqs.  (S8)  in  lieu  of  the  more  costly  direct  inversion 
methods.  The  scheme  can  be  represented  as  follows: 


(«)C‘  -  ^ 


*f.,i 


wherr  j  is  the  iteration  index.  To  test  the  iteration  procedure,  the  volume  source  driven  analytical 
beuchmark  with  constant  properties 


^x,0)  -  0,  1/2,0)  -  0 


was  used  with  the  solution 


■  J  «»  (B.  X) 


where 
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B«-(2a+l))r  . 

Also  of  note  is  the  approach  to  the  "steady  state”  or  equilibrium  solution 

^x)  -  1  -  cosh  (x)/cosh  (1/2)  (63) 

as  s  approaches  mfinity. 

Table  3  displays  the  relative  error  for  increasing  iteration  index  with  opbOJ  at  several  points 
within  the  slab  and  for  a  uniformly  distributed  source  emitting  electrons  at  all  path  lengths.  A 
rather  large  number  of  iterations  are  required  to  achieve  acceptable  accuracy  which  is  considered  to 
be  10**’  for  the  electron  studies  presented  here.  The  approach  to  the  equilibrium  solution  is  shown 
in  Figs.  6a«b  with  a  graphical  comparison  of  the  exact  and  finite  difference  (FD)  solution  given  in 
Fig.  6c  for  As  •  0.01  and  Ax  <■  O.OS.  Again  acceptable  agreement  is  noted.  Table  4  gives  the 
relative  error  variation  with  As  at  large  values  of  s  as  ^s)  approaches  the  equilibrium  solution  for 
Ax  «  O.OS.  There  is  some  indication  that  refining  As  without  refining  Ax  does  not  improve  the 
results.  This  is  an  indication  of  a  limitation  resulting  from  the  iteration  procedure.  Figures  7a,b 
demonstrate  that  in  general  the  best  choice  of  a  is  This  choice  of  a  apparently  combines  the 
accuracy  of  the  explicit  method  and  the  stability  of  the  implicit  formulation.  A  comparison  of 
results  for  different  As's  is  shown  in  Figs.  8a,b.  The  accuracy  is  apparently  limited  by  Ax  as  was 
similarly  found  in  analyzing  the  approach  to  the  equilibrium  solution  (see  Table  4).  Figure  9  shows 
the  improvement  as  Ax  is  reduced  for  a  fixed  As.  Again  the  results  are  not  expected  to  improve  as 
Ax  is  further  reduced  unless  As  is  also  reduced, 
b.  Direct  Matrix  Inversion 

The  iteration  process  can  be  avoided  altogether  by  utilizing  a  tridiagonal  solver  (TDS)  to 
perform  the  required  matrix  inversion  of  Eq.  (S9)  in  the  form 
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X  .  a  A®  -  (I-a)  A-T"  +  Q“  (64) 

That  the  iteration  is  responsible  for  an  apparent  instability  when  Ax  is  reduced  with  As 
fixed  is  clearly  demonstrated  in  Figs.  10a,b.  In  this  figure  and  those  to  follow,  is  the  number  of 
interval  partitions  of  the  slab  length.  The  iterative  solution  oscillates  widely  about  the  exact  solution 
for  small  values  of  As;  whereas,  the  TDS  solution  does  not  and  consistently  provides  accurate  results 
for  all  values  of  Ax.  Thus  there  is  a  compelling  argument  to  abandon  the  iterative  procedure  in 
favor  of  the  TDS. 

For  the  comparisons  to  follow,  a  global  figure  of  merit  to  measure  the  error, 

itb  Ly  ^ 

Egl(*)  •  J  ^  ^  AXi|^(s)  -  ^(s)|  ,  (65) 

i-1 

has  been  found  to  give  a  reliable  characterization  of  solution.  is  the  total  number  of  spatial 
intervals  between  a,b  and  represents  the  exact  solution.  The  variation  of  Eqi,  with  s  for  several 
is  displayed  in  Figs.  lla,b  again  effectively  indicating  the  increased  stability  of  the  tridiagonal 
matrix  inversion. 

To  further  demonstrate  the  superiority  of  the  TDS,  several  additional  test  problems  with 
analytical  solutions  have  been  considered.  The  test  cases  have  been  designed  to  span  the  types  of 
situations  anticipated  to  occur  when  evaluating  electron  penetration  in  actual  electronic  devices. 

The  second  problem  considered  is  for  a  slab  with  constant  properties  p  and  g  and  is  driven 
by  a  boundary  source  The  analytical  solution  is 
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» (h 


-  nx/a,  Smq/p  . 


(66) 


Figures  12  again  show  unphysical  behavior  of  the  iterative  solution  as  Ax  is  reduced.  The  more 
well  behaved  solution  obtained  by  the  tridiagonal  solver  is  shown  in  Fig.  12b.  The  density  near  the 
source  boundary  seems  to  become  less  accurate  as  Ax  decreases.  This  inaccuracy  can  be  eliminated 
if  a  finer  mesh  is  specified  near  the  boundary  than  for  the  interior.  Figure  12c  shows  the 
inprovement  when  only  a  single  small  additional  interval  is  introduced  at  the  boundaries — the 
discrepancy  is  completely  eliminated.  Hgures  13a,b,c  give  the  corresponding  global  errors  for  the 
second  benchmark.  The  global  errors  for  the  variation  of  As  in  the  iteration  and  matrix  solutions 
are  shown  in  Fig  '14a,b.  Again  the  improvement  ^provided  by  the  tridiagonal  solver  is  clearly 
evident 

As  a  final  demonstration,  the  diffusion  coefficient  p  in  problem  1  was  varied  from  0.1  to  SO 
and  the  global  error  for  the  TDS  is  presented  in  Fig.  IS.  The  accuracy  decreases  as  p  increases  most 
likely  resulting  from  the  increased  emphasis  of  the  discrete  approximations  used  for  the  spatial 
derivatives  as  p  increases.  These  results  indicate  that  variable  electron  properties  can  easily  be 
accommodated  with  acceptable  accuracy. 

S.  Dose  Calculation 

In  an  attempt  to  apply  the  discretized  solution  to  a  realistic  problem,  the  dose  in  thick 
aluminum  was  determined.  The  diffusion  equation  describing  the  collided  density  f^  is 

[s  - '  °  fc(«) 


(67a) 
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fe(x,0)  -  0  ,  fe(r^)  -  0  (67b) 

where  r,  is  the  electron  range  in  aluminum  and  all  properties  are  assiuned  constant  with 

D  -  A*/3  . 

The  uncolUded  electron  density  for  normal  incidence  is  given  by 

fo(Jt,s)  -  m  . »?  ■  x/s  .  (68) 

The  electron  density  is  therefore 

^X4)  -  fo(x^)  fe(x4)  (69) 

yielding  the  dose 


for  a  constant  stopping  power  and  the  dose  normalized  at  x  >  0,  the  results  are  shown  in  comparison 
with  other  theories  in  Fig.  16.  Surprisingly  good  agreement  is  achieved  which  is  most  probably 
fortuitous  given  the  crude  nature  of  the  dose  approximation.  The  results,  however,  are  encouraging. 


B.  2-D  Diffusion  Theory 

The  diffusion  equation  in  two-dimensions  is 
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^  (x,y4)  -  ^  p(x,y^)  ^  ^x,y^)  + 

+  ^  P(  W)  ^  ^x.y4)  -  q(x,y^)  ^x,y4)  +  S(x,y4) 

with  boundary  and  initial  conditions 

^a,y4)  ■  oi.  .  ^b,y^)  -  (Tr 
^X,C4)  ■  OSj.  ,  ^x4«s)  -  <^3 
^x,y,0)  -  f(x,y)  . 

Integrating  over  cell  (ij)  depicted  in  Fig.  17  and  employing  the  approximations 
dx  dy  ^x,y^)  -  ^  (Ayj  +  Ayj+j)  (Ax,  +  Ay,^i) 

JAy 


(71a) 


(71b) 

(71c) 

(71d) 


f  ha  p  at-  *IJ-*LI-1 

Jab  «  Ay,  [  2  J 

results  in 

d  A  • 

-  laij  +  Cy  +  by  ^]  + 

+  fy  ^J+I  +  «y]  ■  5lj  ^  ^ 


(72) 


where  %j.  by.  c,j,  djj,  ey,  fy,  ?i,j  and  Qj  are  ceil  constants.  Several  differencing  schemes  for 
the  path  length  variable  will  now  be  discussed:  Note  that  one  could,  in  addition,  apply  the 
continuous  analytical  continuation  techniciue  (§in.A.2),  however,  that  approach  win  only  be  applied 
in  the  3-D  study. 
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1.  Two-Step  Alternating  Direction  Implicit  (ADI) 

In  this  scheoM,  the  path  length  derivative  is  split  in  the  interval  As  and  the  resultant  set  of 
equations  is  first  solved  in  the  x  direction  and  then  in  the  y  direction.  In  step  1, 


ds  ”  As/2 

yielding 

+  Cy  j  + 

(73) 

+  ^4  ^  ^  ^-1  ^+1  ®u  " 


As 


i 

3 


Eq.  (73)  is  solved  using  the  TDS  in  the  i  direction.  The  second  step  consists  of  the  approximation 


(74) 


where  the  terms  in  the  brackets  are  the  same  as  in  Eq.  (73)  evaluated  at  die  indicated  discrete  path 
length.  The  TDS  is  then  applied  in  the  j  direction  giving  the  final  solution  at 
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2.  Iterative  Mediods 

As  io  the  1-D  case,  tiiere  exists  an  iterative  formulation  of  the  solution  to  Eq.  (72)  when  the 
backward  finite  difference  approximation  to  the  path  length  derivative  is  employed.  The 
formulation  becomes 


<r‘  -  «2|.  <i«l  (73» 

with 

^  (5)  -  a  ♦  (1-a)  4{|i.  0  i  a  i  1  .  (75b) 

In  the  following,  several  of  the  many  possible  ways  of  solving  for  are  indicated: 

a.  Line- Jacobi  (U)  with  TDS 
.  (D^ 

P«) 


b.  Succasstv,  OvwralixatioD  <SOR)  with  Itontion 

-  «  Fl«4Ki,  w^.\,  (k'DCij.  +  (!-")  ™ 

0  <  «  <  1  (77b) 

c.  Successive  Line  Overrelaxation  (SLOR)  with  TDS 

(78a) 

(k)^»  .  u  ♦  (l-«)  (78b) 

0  <  w  <  1  .  (78c) 


3.  Test  Results  and  Discussion 

A  source  driven  test  problem  with  constant  properties  and  a  uniform  source  of  the  form 
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s(s)  -  Qoe-^ 

will  be  wed  to  test  the  various  solution  procedures  described  above.  For  zero  density  on  the 
boundary,  the  solution  is 

^x,y,s)  -  ^ 

m,n 

with 

cos  (B„  x)  cos  (B^  y) 

Bb*  ■  (2n+I)  T/a 
B„^  -  (2m+l)  T/b  . 

A  more  numerically  amenable  solution  is  obtained  by  rearranging  the  series  in  £q.  (79a)  in  the  form 

oo  n 

0(x,y4)  -  ^  ^  Ab,^.«  .  (80) 

naO  maO 

Figures  18a,b,c  show  comparisons  of  the  ADI,  LJ  and  SOR  solution  at  y  >  0  and  s  >  0.5  to 
the  benchmark  solution  for  the  variation  of  Ax  and  Ay.  ^  ^  number  of  intervals  in 

the  X  and  y  directions  respectively.  All  methods  seem  to  give  comparable  accuracy.  The  ADI  and 
SOR  methods  require  the  largest  computational  time  with  the  L)  technique  requiring  the  least 
Further  testing  will  determine  the  final  selection  of  the  best  method  for  the  anticipated  vulnerability 
stixUes. 

4.  Transport/DilTusion  Coupling 

To  demonstrate  the  coupling  of  a  transport  and  diffusion  calculation.  Fig.  19  shows  a 
simuhited  two  dimensional  approximation  to  a  one-dimensional  transport  result  For  this  calculation. 


(79a) 


(79b) 
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the  density  profile  at  s  •  5  was  taken  from  an  analytical  transport  theory  calculation**  and  used  as 
the  initial  condition  to  the  2-D  LJ  diffusion  solution.  Both  the  analytical  and  diffusion  theory 
«iigniariftn<t  were  then  run  to  s  15  and  Fig.  19  shows  the  comparison  for  several  Ly. 
Acceptable  agreement  is  seen  except  possibly  near  the  peak.  The  comparison  can  be  improved, 
however,  by  taking  a  finer  mesh  in  region  of  the  peak.  Additional  studies  have  been  performed 
inHiraritig  that  the  Coupling  of  transport  theory  and  diffusion  theory  in  the  above  manner  is  a  viable 
ffaiffiiiatinnai  method.  It  is  thk  coupling  that  will  serve  as  the  basis  for  future  development  of  the 
diffusion  theory  model. 


IV.  Development  of  analytical  Benchmark  Methods 

The  trademark  of  the  electron  transport  theory  methods  devetopment  at  the  University  of 
Arizona  has  been  and  will  continue  to  be  the  pursuit  of  analytical  benchmarks  while  in  the  process 
of  forming  numerical  algorithms.**"**  The  interaction  of  these  two  components  has  proved  to  be 
invaluable  in  the  development  of  reliable  and  efficient  numerical  methods. 

The  benchmark  ^methods  to  be  presented  below  represent  an  entirely  new  concept  in 
benchmarking.  While  the  benchmarks  established  thus  far  are  not  directly  applicable  to  the  electron 
transport  application,  the  methods  will  certainly  prove  helpful  in  developing  new  benchmarks  that 
wiU  be. 

A.  Benchnuurks  from  Laplace  Transform  Inversion 

The  method  in  this  section  is  taken  from  Ref.  (30).  Recently  a  highly  reliable  numerical 
Laplace  transform  inversion  has  been  developed  based  on  the  transformation  of  the  Bromwich 
inversion  integral  into  a  cosine  integraL  The  cosine  integral  is  then  represented  as  an  infinite  series 
of  integrals.  These  integrals  are  evaluated  using  a  Romberg  integration  scheme  and  the  convergence 
of  the  series  is  Mcelerated  with  an  Euler-Knopp  acceleration  algorithm.  In  addition,  the  inversion 
is  performed  on  several  contoun  and  agreement  is  required.  A  demonstration  of  the  accuracy  of 
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the  inversion  is  shown  in  Fig.  20  where  the  algorithm  accuracy  c  was  varied  from  10~*  to  10~*. 
This  algorithm  will  now  be  used  to  provide  highly  accurate  (3-S  digit)  benchmarks  for  the  stationary 
one  velocity  Boltzmann  equation. 

Consider  the  one-group  Boltzmann  equation  in  a  semi-infinite  medium  with  isotropic 
scattering 

dM'«x,M')  .  (81) 

1 


for  n  <  0.  By  letting  -  1/s  where  s  is  a  complex  variable,  we  note  that  Eq.  (82)  becomes  a 
Laplace  transform.  Thus  upon  inversion 


The  inversion  can  be  performed  numerically  since  ^0,-m)  is  known  to  be 


^0,-m)  -  1  -  H(m) 


(84) 


-  48  - 


The  analytical  continuation  of  ^0,-m)  to  the  complex  s  plane  can  easily  be  done  thus  allowing  the 
numerical  inversion. 

The  numerical  inversion  is  applicable  to  any  linear  transport  problem  in  which  the  exiting 
distribution  is  known.  Figure  21a-d  show  results  for  several  basic  transport  problems  with  isotropic 
scattering  in  a  semi-infinite  medium  including 

1.  beam  source  illuminating  the  boundary 

2.  one-half  space  Milne  problem 

3.  uniform  source  in  right  half  space. 

The  method  is  equally  applicable  to  the  case  of  anisotropic  scattering. 

In  the  future  an  attempt  will  be  made  to  apply  the  numerical  inversion  to  the  discretized 
2-D  transport  equation  in  order  to  treat  the  s  variable  continuously  and  provide  a  reasonable  2-D 
transport  benchmark  for  electrons. 

B.  Moments  by  Continuous  Analytical  Continuation 

Accurate  moments  of  the  Spencer-Lewis  equation  in  one-dimension 

^  +  M  ^  +  A(s)j  <Kx,tiA)  -  A(s)  I  dn'  fOt»  «xX^)  (86) 


in  infinite  geometry  are  of  considerable  importance  in  the  reconstruction  of  the  density.  The 
moments  are  defined  by 
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M^(s)  a  exp 

where  is  the  Legendre  polynomial  of  order  L  satisfies  the  following  set  of  ordinary 
differential  equations 

H  +  hi  -  Hs)  (l-u,i)  +  Q^(s)  (88) 

where 

a  I  dp'  P^(p')  f(p',p) 


(S  .1  foo 

ds'  A(s')  dp  ?iOi)  dx  x“  ^x,p^) 

0  J-1  J-oo 


(87) 


and  contains  boundary  densities  known  from  an  independent  calculation  or  approximated  from 
a  lower  order  theory.  Then  by  expanding  Q(^  and  in  Taylor  series  in  the  interval  s^.^  <  s  < 

oo  y 

n>4) 

n*0 

oo 

M4k(s)-  ^  ^(s-Vi)“  . 
n»0 

we  find  the  recurrence  relation 


■  -  a£  -  hi  + 


(89) 
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j-0 

MflAk  “  • 

M/^(s)  can  then  be  accurately  obtained  by  defining  as  many  s  intervals  as  required  for  the  desired 
convergence. 

The  density  can  be  reconstructed  in  the  usual  way  by  an  expansion  in  orthonormal 
polynomials 


(90a) 


(90b) 


(90c) 


One  must  be  alert  to  the  possibility  of  round  off  error  contamination  in  Eqs.  (89)  and  (90c). 
Equations  (90)  are  expected  to  provide  accurate  results  away  from  the  s,x  wavefront.  Acceleration 
techniques  are  then  applied  to  Eq.  (90a)  in  an  attempt  to  mitigate  round  off  error  and  provide  more 
accurate  solutions  near  the  wave  front  discontinuity. 
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Table  2  Comparison  of  results  obtained  for  particle  streaming  through  a  vacuum  from 
the  (SRI  method  (top),  the  (SR)  method  with  ray  mitigation  (middle),  and  the  (SCSN) 
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Table  3 


Iteration 

Error 

s  «  0.01 

J 

*4 

*e 

2 

1.0000E-tO0 

1.0000E400 

3 

4.4420E-01 

4.4420E-01 

4 

2.6199E-01 

2.6199E-01 

5 

1.5484E-01 

1.7313E-01 

6 

1.1012E-01 
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Figure  1.  Angular  cone  construction  for  i  >  t'  + 1.  Cross  hatched  area  is 
the  extended  region  of  influence  for  direction  m.  Double  cross  hatched  area 
is  that  fraction  within  the  target  cell  ij. 


Figure  2.  Angular  cone  construction  for  »  <  i*  +  l.  Cross  hatched  area  is 
the  extended  region  of  influence  for  direction  m.  Double  cross  hatched  area 
is  that  fraction  within  the  target  cell  ij. 


ScMlar  Flux 


solution.  Obuined  firom  th®  coarentionnl  Sn  (DD)  method  t 
sMthod  for  Bne  (3.7  mfy)  ia«b  ceU*.  a)  (DD)  JW  *  4,  b)  (SCSNl  Af  - 
c)  (DD)  Af  *  8.  d)  (SCSN)  M  *  8,  ®)  (DD)  M  *  12,  f)  (SCSN)  M  "  12.  ^  “ 
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Table  4 

Approach  to  "Steady  State*  ^x)  *  0.11318 


As 

_s 

^04) 

ERROR 

0.01 

OJ 

0.11229 

7.9  X  10-1 

0.02 

1.0 

0.11277 

3.62  X  10-» 

0.04 

2.0 

0.10913 

3.58  X  10-s 

0.08 

4.0 

0.09972 

1.19  X  10-1 

0.02 
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Fig.  6c  Comparison  with  EXACT  solution 


Fig.  7a  Comparison  for  different  o's  (As  =  0.01,  Ax  =  0.03)  for  (a)  x  =  0  and  (b) 
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10l»  Variation  of  Ax:  Trldlagonal  —  source  driven  benchmark 
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Fig. 


Fig*  I'ia  Global  error  cooiparlson  with  Prob  2 


Fig.  13b  Global  error  comparison  with  Prob  2:  Trldlagonal 


Global  error;  Trldlaaonal  -  Prob  2 
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FJg.  l<ib  Variation  of  As 


Fig.  16  Dose  comparison  In  aluminum 
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I'ig*  18a  Comparison  with  2-D  benchmark:  Ax,  Ay  senslcivlcy-AUl 
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Fig.  20  Variation  of 


(uo  -  1) 


